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1. INTRODUCTION 


Advanced composite materials, especially graphite/epoxy, are being 
applied to airplane structures in order to improve performance and save 
weight. An important consideration in composite design is the residual 
strength of a structure containing holes, delaminations or interlaminar damage 
when subjected to compressive loads. While elastic behavior of composites has 
been studied extensively in recent years, the viscoelastic response of these 
materials is not well understood. Recent studies by several investigators 
have revealed the importance of viscoelastic effects in polymer-based 
composites [1-4]. The viscoelastic effect is particularly significant at 
elevated temperature/moisture conditions since the matrix material is strongly 
affected by the environment. 

The solution of viscoelastic problems in composites has been limited to 
special cases which can be solved by classical lamination theory [2-4], In 
this report, a finite element procedure is presented for calculating 
time-dependent stresses and strains in composite structures with general 
configurations and complicated boundary conditions. Using this procedure the 
in-plane and interlaminar stress distributions and histories in notched and 
unnotched composites have been obtained for mechanical and thermal loads. 

Both two-dimensional and three-dimensional viscoelastic problems are analyzed. 
The effects of layup orientation and load spectrum on creep response and 
stress relaxation have also been studied. 
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2. CONSTITUTIVE EQUATIONS FOR ANISOTROPIC VISCOELASTIC MATERIALS 

Consider a linear anisotropic viscoelastic material subjected to both 
mechanical and thermal loads. Using the contracted notation for the stresses 
and strains, the constitutive relations may be written in the following 
integral form [5]: 

t SI,(T) 

o,(t) =/ C. ,(T.t-T) dt (1) 

a. CO 

with *£ =£.- £.* 

J J J 

In the above, cl and are the stress and the total strain components along 
the principal material directions (123 axes). C, . are the rel axation modul i , 

* J 

T is temperature and t is time, e.* is the thermal strain component in the 
stress-free state and can be expressed as a function of the coefficient of 
thermal expansion aj and the stress-free temperature T* by 

e j* ■ £ a j< T » dT w 

In this study the aj is assumed to be independent of time and temperature, 

that is, e.*= a. AT with AT = T - T*. 

J J 

The basic viscoelastic properties C. . can be determined from the 

J 

experimental characterization of a unidirectional composite at various 
temperatures over a long period of time. However, this technique is quite 
time consuming and an alternate technique based on the time/temperature 
superposition is frequently used to minimize the test time required for 
polymer based composites. The alternate technique consists of testing a 
series of specimens at higher temperatures. The test results are plotted and 
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the resulting curves are then shifted horizontally (and sometimes also 
vertically) with respect to a reference temperature to obtain the master 
relaxation curves for the C. . components. The amount the curves are shifted 

1 J 

is termed the shift factor <J>-j j . Materials whose viscoelastic properties can 
be characterized using this procedure are referred to as thermorheol ogical ly 
simple materials. For this special class of materials the relaxation 
functions C.. in the principal material directions may be represented in the 

* J 

following form 

Ci j(T ,t) = C-j j(T 0 , j ) (3) 

where T 0 is the reference temperature for the master curves and 5-j j is the 

reduced time. The reduced time is related to the shift factor <Hi(T) by 

t 

? ij(t) = / *ij (T(s) ) ds (4) 

o 

As in the elastic case the relaxation functions TT-jj along the load 
directions xyz (Fig. 1) can be obtained from the C-j j of the principal material 
directions (123) by the use of tensor transformations. With the transformed TT-j j , 
the viscoelastic constitutive equation for an anisotropic material along the xyz 
coordinates becomes 

V*' - -57^ dT (5) 

T 

where ^(t) = / ♦- .(T(s)) ds (6) 

1 j o J 
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3. VARIATIONAL THEOREM FOR LINEAR THERMO-VISCOELASTICITY 


A variational theorem for linear viscoelastic materials is given by 
Christensen [6] and is extended here to include noni sothermal effects. For an 
anisotropic material the functional * can be defined as 


s=t 

» - / / 

V s=-‘ 


T=t-S 


/ {j C- .(T,t-S-x) 

T— — OO V 


( T ) 


X= - <» 


8 e (s) . 3 e (s) 

■ CjjfT.t-S-t) yj- } dTdsdV 


- / / ' 1 T (t-s) |Hi (s) ds dA 

S(j S = -“ 


(7) 


where V is the volume domain, u-j is the displacement, and S a is the portion of 
boundary where the traction T-j is prescribed. 

Taking the first variation of n and letting Sn = 0 along with the 
commutative relations of Stieltjes convolution yields 


6tt 


J S_t ( / Oj(t-s) iiaisl dV - / Ti(t-s) dA ) ds = 0 

s=-“ V s S a s 


( 8 ) 


The solution to the quasi-static thermo-viscoelastic boundary value problem is 
found by extremizing the functional ir [6,7], This variational method will be 
used in the following finite element formulation. 
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4. FINITE ELEMENT FORMULATION 


Consider a symmetric laminate with arbitrary ply orientations subjected 
to both mechanical and thermal loads. It is assumed that the prescribed 
temperature is uniform throughout the laminate at any time instant and that 
the applied surface traction is separable function of position and time. 

These assumptions allow the displacement solutions to be expressed as the 
product of two separate functions, one involving only spatial coordinates and 
the other involving only spatial coordinates and the other involving the time 
variabl e. 


In the finite element development, 4-node quadrilateral 2-D elements and 
8-node 3-D solid elements are considered. Using an isoparametric formulation 
the displacement fields within an element are interpolated as 

{ u } = [ N ] { q } (9) 

where { q } is a nodal displacement vector whose components are a function of 
time only. The strains { e } can be obtained by the differentiation of the 
displacement field in Eq. (9) and are expressed as 

{ * } - C B ] { q } (10) 


Use of Eqs. (9) and (10) and the variational theorem discussed in Eq. (8) 
leads to the following equilibrium equation for an element 


L Vli"V T ’ l - T) B Jn t “ ! ^ a<T) dt 


k B im < T - ‘-0 Tf <T)d ' d(1 *f T „ \m T k(t) dr 


(ID 


where fl is the volume of the element. 
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Note that the first term on the right hand side of Eq. ( 11 ) is the force 
vector due to thermal loading and the second term is the reactive force 
vector. When Eq. ( 11 ) is assembled the reactive force vector vanishes 
everywhere except on the boundary where the traction is prescribed. 

Before proceeding further, the following abbreviations for C-jj of a 
unidirectional ply are defined as 


C 1 = Cn, C2 = C12, C3 = C13, C4 = C22. C5 = C23, 
C6 = C33, C7 = C44, Cg = C55, Cg = Cgs 


( 12 ) 


These nine relaxation functions are expanded in terms of exponential 
series so that the integral equation ( 11 ) may be easily calculated [ 3 , 7 ]. The 
exponential series is 


NT 



C i (t ) = C-j 0 + I ^io) ex P 

( — t / X-j (j) 

( 13 ) 

(0=1 

- 


where Aj^ is the relaxation time obtained 

from the master curve. The TTfj 

i with 


respect to the xyz axes can be obtained with the use of tensor transformation. 
The results are 

9 

C i j ( t ) = I n i jr Cp(t) (14) 

r=l 


where j r is the component of the transformation coefficient matrix for each 
layer. Note that n jj r is symmetric with respect to indices i and j. The full 
expressions for n ij r are given in Appendix A. 

From Eq. ( 11 ) the equation for the element stiffness matrix k mn of an 
element is 
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kmn U" T ) = !q B im Cjj (T, t - t) Bj n dQ 


y nt 

= I { kmnr.o + Z ^mnr,w ex P E-(?r “ £r)/V«3 } (15) 

r=l u=l 

where k mnr>a) is 

kmnr,w = C ra > J fi Bj m nij r Bj n dfl (16) 

Similar to the elastic case, the element stiffness matrix and force vector are 
assembled over the whole domain to yield the following global equation for the 
displacement u n (t): 

Kmn ( 0 u n (0) + / K mn (? - C") — — P-( T ) dT = F m (t) + F m (t) (17) 

0 3t 

where 

9 NT 

K mn(s ■ S') = I { Kmnr.o + I, Kmnr,o> exp[-(c r - s r )/Vu>] } (18) 

r=l (d=1 

and m, n = 1, 2, 3, , NDT, NDT = total number of degrees of freedom. 

In the above, K mnrjW is associated with the global stiffness matrix, 
F m t (t) is the component of the global residual force vector due to thermal 
load and F m r (t) is the component of the global force vector due to the 
prescribed tractions. 

A direct integration of Eq. (17) requires enormous computer storage space 
for the stiffness matrices and the displacement vectors of previous times and 
is not feasible in most problems of practical interest. To overcome the 
storage limitations and reduce computing time, a numerical scheme similar to 
that employed by Taylor et al . [7] for an isotropic material is used. Eq. 

(17) can then be replaced by a summation of the integrations over a series of 
time intervals At. Within each subinterval of time the dependent variables 
u n (t) may be approximated by a linear variation. That is. 
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> (tl 5 fortj.j <t <tj (19) 

where Au n (tj) = u n (tj) - u n (tj_ 1 ). (20) 


Using the above approximation for the derivatives of u n (t), Eq. (17) at 
t = tp becomes 


p t i 9 NT , x 

If I {^mnr,o + I ex P C - (?r>p “ 4r)/VoJ} dt — Or- . J- 

j=l tj_i r=l oi=l J 

( 21 ) 

9 NT 

= ^m^(tp) + f r m r (' t p) ~ I {^mnr s o + I ^mnr.to exp(-c r ,p/Vw) } u n(0) 

r=l tu=l 

where 4 fjP = ^ r (t p ). 


Finally Eq. (21) leads to a set of . algebraic equations 


9 NT 

I ( K mnr,o + l K mnr,u h raJ (At p )} Au n (t p ) = F m t (t p ) + F m r (t p ) 
r=l w=l 

( 22 ) 

9 9 NT 

I K mnrj0 u n(Tp-l) - 1 l 9mr,u>(t p ) 

r=l r=l u=l 

where 

i t j ^ 

hp^lAtj) = -j ^-7 / sxp[-(t r » > j - ? r)/Va>3 dT (23) 

J \j-l 


9mr,w(T p ) = Kmnr.u ex P( - ?r,p/ Va>) u n(0) 

(24) 

P-l 

+ I ^mnr,£D ex PC"(V,p ~ V ,j )/ Vw3 h ra ,(Atj)Au n (tj) 
j = l 


If the temperature is constant within the time interval Atj, h rw (Atj) can 
be evaluated exactly as follows 
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hp<jo( Atj ) - ^poj[l ~ 6Xp( -A^p ^ j/Xpm) ] /A^p^j 


(25) 


where A? r,j = 4 r ,j “ 4 r ,j-l (26) 

If the temperature is not constant during Atj, the shift factor 
corresponding to the specific temperature variation within the time interval 
can be used to obtain the reduced time from Eq. (4), from which h ra) (Atj) can 
be determined numerically. 

Using Eq. (24), a general recursive expression can be derived for p > 1. 

The recursive formula is 

9mr,o)(t p ) = exp(-Ac r>p / V<o) E9mr ,ui(tp-l ) + K mnr,o» h ra) ( Atp_l.) Au n(tp-l)J (27) 
with 

9mr , a)(to) = °> h rJ A t 0 ) = 1, and Au n (t 0 ) = u n (0). 

The residual thermal force vector can be updated using a numerical 
algorithm similar to that used for the stiffness matrix. The element residual 
thermal force vector represented by the first term on the right-hand side of Eq. 
(11) is expressed as 

9 NT 

fm^(t) = f Jq ^im I ^ijr l^ro + I C ru) exp[-(? r - zf)/\ r oJ} — ~ dft dx 

r=l oj=1 J 

(28) 
(29) 

The oj are the coefficients of thermal expansion along the xyz axes. 

These force vectors can be assembled to obtain the global residual thermal 
force vector F m t (t) in Eq. (21). 


9 NT 

= I 7 {fmro + I fmru ex p[ _ (?r ~ — fa “ dx 

r=l - 00 u=l 


fmrco = Bim (^ijr ^ro) a j) dn 
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For the non-i sothermal case, the solution procedure is to descretize 
temperature histories so that the temperature change occurs only in elastic 
time step A-A 1 (At = 0) with the temperature remaining constant during the 
subsequent viscoelastic time step A' -B (See Fig. 2). Within the time step 
A-A 1 , the reduced time is not changed and h rw (Atj) = 1. During the 
viscoelastic time step A'-B, the difference in the reduced time A? r (tj) is 
determined from the time difference Atj multiplied by the corresponding shift 
factor. Eq. (25) is then used to calculate h ra) (Atj). 

After Eq. (22) is solved for Au n , the nodal displacements at current time 
t p are obtained from the relation u n (t p ) = u n (t p _i) + Au n (t p ). Once the 
displacements are found, the strain field can be determined. Finally, the 
viscoelastic constitutive relations in Eq. (5) are used to calculate the 
stresses. 

It is noted that the above formulation is valid for general 3-D 
viscoelastic problems. In the special case of a 2-D analysis, the "C-jj matrix 
of each layer is reduced toTJ-jj and is averaged through the laminate thickness 
to obtain the extensional stiffness matrix A-jj. Also note that the 
approximation of the dependent variables u n (t) by a linear Lagrangian 
interpolation function may cause significant error accumulations since the 
solution at the current time is affected by the previous solutions. Such 
error is expected to grow as the number of solution iterations increases. In 
order to minimize the error accumulation, higher order interpolation functions 
are needed. However, such a formulation becomes extremely complex and 
cumbersome since more than one set of previous solutions are required to 
obtain the current solutions. 
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5. NUMERICAL RESULTS FOR THE TIME-DEPENDENT RESPONSE 


Numerical results on time-dependent stress/strain fields are obtained for 
graphite/epoxy composites using the formulation derived in the previous 
sections. Elastic material properties used in the analysis are those typical 
of graphite/epoxy composites as shown in Table 1. With these engineering 
material constants, the elastic stiffness matrix [C-jj(O)] can be obtained. 
Viscoelastic relaxation functions C-jj(t) are then generated by multiplying the 
elastic stiffness C j ( 0 ) by specific time varying functions. In this study, 
the C n is assumed to be independent of time and temperature while other 
relaxation functions are assumed to have the same time-varying function f(t). 
The function f(t) is taken from Flaggs and Crossman's experimental curve 
[3,8] and is expressed in terms of an exponential series containing 11 terms 
(See Table 2). The shift factors corresponding to various temperatures as 
given in [8] are the same for all Cij(t) master curves. The stress free 
temperature for the graphite/epoxy laminate is assumed to be 350°F. The shift 
factors at various temperatures are shown in Table 3. 

5.1 Unnotched Laminates 

The first analysis is to verify the accuracy of the present approach by 
comparing solutions with the incremental classical lamination theory (CLT) 

[2]. The case studied is a ( 0/45/90/-45 ) s quasi-isotropic laminate (1 in. 
wide by 2 in. long) subjected to thermal loads. The mesh pattern used in the 
finite element solution consists of 50 4-node square elements with a total of 
132 degrees of freedom (Fig. 3a). The compliance function S 22 (t) and S g6 (t) 
at T = 75°F are given below: 

S 22 = 0.7143 x 10" 6 + 0.00385 x lO" 6 t 0 * 33 (psi ) _1 

S 66 = 1.1111 x 10- 6 + 0.00680 x IQ" 6 tO- 31 (psi)" 1 
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where t is in minutes. Other constants are assumed to be independent of time, 
i.e. E u = 18 x 10 6 psi and v i2 = 0.34. The normalized master curves Q 22 (t) 
represented by the exponential series and the power law are illustrated in 
Fig. 4. 


Using these material properties the stress and strain histories have 
been calculated for two loading conditions. Case I involves an instantaneous 
change of temperature from 350°F to 75°F at time t = 0 with the temperature 
held at 75°F thereafter. Case II is a sudden temperature change from 350°F to 
160°F at t = 0 with the temperature kept at 160°F afterwards. The in-plane 
strains e x obtained for Case I are shown in Fig. 5. It is seen that the 
present finite element solution agrees well with the incremental solution from 
the CLT approach, the maximum discrepancy being 1.5%. Note that the strain e x 
decreases with time. After 1.2096 x 10 6 seconds (two weeks), the magnitude of 
the strain decreases by 10.5% in Case I (AT = -275° F and T = 75°) and 42.7% 
in Case II (AT = -190° F and T = 160°F). The viscoelastic stress o x in the 0° 
layer of the (0/45/90/-45) s laminate, normalized with respect to the initial 
thermal stress at time t = 0, is shown in Fig. 6 for both cases of loading. 

The initial thermal stress o x (0) in the 0° layer is -5380 psi for AT = -275°F 
(Case I) and -3717 psi for AT = -190°F (Case II). It is noted that while the 
elastic solution is linearly proportional to the temperature change, the 
stress in a viscoelastic analysis relaxes with time in different proportions 
(see Fig. 6); at the completion of 336 hours the thermal stress o x relaxes 
11.9% in Case Moading and 48.4% in Case II. The difference in the amount of 
stress relaxed is due to the fact that relaxation moduli TJ-j j depend on both 
the current temperature and time. The relaxation moduli decrease more 
significantly with time at the higher temperature (T = 160°F, Case II) than at 
the lower temperature (T = 75°F, Case I). Comparison of these results also 
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shows the nonlinear effect of temperature on stress relaxation. 

Figure 7 depicts the laminate strain history normalized with respect to 
the initial strain at time t = 0. The initial strain values are 
e x = -4.56 x 10"4 in/in for Case I loading and e x = -3.15 x 10"4 in/in for 
Case II. It is interesting to note that in both cases the magnitude of 
laminate strain e x decreases rather than increases with time. After t = 1.2E6 
sec. (two weeks) the magnitude of laminate strain in Case I (AT = -275°F and 
T 0 = 75°F) decreases 10.5%, while the laminate strain in Case II (AT = -190°F 

and T 0 = 160°F) decreases 42.7%. 

The effects of different temperature spectrums on stress relaxation and 
creep strains for (45/-45) s graphite/epoxy laminates were also investigated. 
Two conceivable temperature histories are shown in Fig. 8. The dotted line 
(Path A) in the figure denotes that the laminate is cooled down slowly in a 
stepwise fashion from T = 350°F at t = 0 to room temperature T = 75°F at t = 

98 days. The temperature variations are composed of a series of elastic steps 
(At =0) and isothermal processes (constant temperature). The solid line 
(Path B) indicates that the laminate is cooled down suddenly from T = 350°F to 
room temperature at t = 0, after which it is subjected to a cyclic temperature 
variation. The residual stresses and strains in the 45° layer associated with 
these two temperature histories are presented in Fig. 9 and Fig. 10, 
respectively. Although the laminate eventually reaches the same temperature 
(75°F) after 98 days, the residual thermal stress and strain are much greater 
in magnitude for the cyclic temperature history (Path B) . The sum of thermal 
stresses induced during all elastic step changes of temperature (At =0) is 
equal to the elastic thermal stress due to the net temperature change from the 
stress-free temperature (350°F) to room temperature (75°F). That is, in the 
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Path A spectrum the sum of stresses T xy shown by OA x + l\ 2 h 3 + A 4 A g + A g A ? + 

A gA g + A 1Q A n + A 12 A 13 is equal to the stress T xy denoted by OB L . In Path B 
loading, the amount of stress t xy represented by 0B 1 - B 2 B 3 + B 4 B 5 - B g B 7 + 
B 0 B a - is equal to OB,. However, the temperatures at which 

the stress relaxation occurs are different. At high temperature the rate of 
stress rel axation is high, hence the laminate undergoes more relaxation in 
Path A spectrum. For example, during the first time interval the stress T xy 
relaxes from Aj to A 2 in Path A spectrum while T xy relaxes a smaller amount 
from B 3 to B ? in the Path B spectrum. Thus, the stress histories depend 
strongly on the specific load spectrum applied. 

Another example considered is the time-dependent response of a symmetric 
laminate at room temperature (75°F) when subjected to mechanical loads. The 
history of mechanical loading is plotted in Fig. 11. The average stress 
applied is equal to 1,728 psi initially and is held constant until t = 24 
hours, after which the laminate is unloaded elastically. The creep and 
recovery behavior for the three different laminates, (0/90) s , ( 45/ -45 ) s , and 
(0/45/90/-45) s layups, are plotted in Fig. 12. As expected, the (45/-45) s 
layup exhibits the most significant creep response among these three 
laminates. The initial in-plane strains e x along the loading direction are 
equal to 0.564 x 10“3 for the (45/-45) s laminate, 0.24 x 10"3 for the 
( 0/45/90/-45 ) s laminate and 0.177 x 10"3 for the (0/90) s laminate. After 24 
hours the strain e x increases by 10.75% in (45/ -45 ) s laminate, 1.16% in the 
(0/45/90/-45) s laminate and 0.5% in the (0/90) s laminate. When the load is 
released in the elastic step (At =0) at t = 24 hours the strain e x is 
suddenly reduced to 0.607 x 10"4 in the (45/-45) s laminate, to 0.279 x 10~5 in 
the (0/45/90/-45) s laminate and to 0.849 x 10“6 in the ( 0/90 ) s laminate. Note 
that the creep recovery rate of the (45/-45) s laminate is also much higher 


14 



than the other two laminates. 


It must be noted that a sufficiently small time interval is usually 
required in order to obtain accurate solutions to Eq. (22). The number of 
time steps needed generally depends on the shape of the applied load spectrum 
and more time steps are required for a complex load history. A criterion used 
for choosing the proper time step size is to compare the creep response of a 
unidirectional laminate due to a unit step stress with the time variation of 
the components in the compliance matrix. For instance, if a unit step stress 
0 X is applied the resulting creep strain e x (t) must 136 equal to the compliance 
component Ty^t). In all of the above analyses, 12 time steps with variable 
interval At were used. Typically the At value is set to 100 sec. at t = 0 and 
increases with time to a maximum of 4.3 x 10^ seconds. 

5.2 Notched Composites 

The geometry considered for a notched composite is a 1 in. wide by 2 in. 
long laminate with a circular hole of diameter 0.25 in. at the center. Since 
there are no viscoelastic solutions available for comparison with this case, 
the accuracy of the finite element solution is estimated based on the elastic 
results at time t = 0. To compare the elastic solutions, three finite element 
mesh geometries were used. The resulting elastic stress concentration factors 
are given in Table 4. All three solutions compare well with the solution by 
Nuismer and Whitney [10], indicating the adequacy of the mesh geometry. As a 
result (of this congruity) an intermediate mesh pattern with 670 degrees of 
freedom is used for the following viscoelastic analysis (see Fig. 3b for 
mesh geometry). 

A uniform stress of a x = 20,000 psi is applied at the remote boundary of 
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the laminate at a constant temperature of 122° F. Any residual thermal stress 
which may exist in the laminate at the completion of the curing process is 
neglected. The in-plane circumferential strains at the hole edge are 
plotted in Fig. 13 as a function of $ for the (45/0/-45/90) s laminate over a 
time period of 10® seconds. Similar results are shown in Fig. 14 for a 
( 45/ -45 ) s laminate. In both cases, the magnitude of the circumferential 

strains e<f, increases with time. At t = 0, the circumferential strain at <f> 

= 90° is 0.897 x 10“® for the (0/45/90/-45 ) s laminate and 0.155 x 10“ 1 for the 
(45/-45) s laminate. As expected the strain in the (45/-45) layup increases 
at a much higher rate (51.3%) than the strain in the quasi-isotropic 
1 ami nate. 

The relaxation of stresses in graphite/epoxy laminates subjected to a 
uniform strain e x = 0.003 in. /in. is also studied. The stress averaged 
through the thickness is obtained as a function of time. The average 
circumferential stress is shown in Fig. 15 for a ( 45/0/-45/90 ) s laminate 
and in Fig. 16 for a (45/-45) s laminate. In the (45/0/-45/90) s laminate, the 
maximum circumferential stress occurs at <t> = 90°. This stress decreases 
from an initial value of 65328 psi to 63584 psi at t = 10® sec. to 61205 psi 
after 10® seconds. In the (45/-45) s layup, the stress distribution is 
somewhat different from the quasi-isotropic case as can be seen in Fig. 16. 

For the (45/-45) s laminate the maximum circumferential stress occurs at about 
<t> = 60° rather than <j> = 90°. At <p - 60°, the stress relaxes from a maximum 
value of 27271 psi to 24102 psi after 10® seconds and eventually to 19472 psi 
after 10® seconds. Additionally, the location of the maximum stress 
concentration moves slightly from <j> = 60° as time elapses. As before, the 
(45/-45) s laminate exhibits a much higher stress relaxation rate than the 
quasi-isotropic layup. 
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5.3 Three-Dimensional Viscoelastic Results 


A. Verification Studies 

Since there is no 3-D viscoelastic solution available for a laminate with 
a circular hole, the two verification studies are limited to elastic case. In 

the first study, a 2 in. by 2 in. by 0.1 in. thick isotropic plate with a 

circular hole of 0.25 in. in diameter is analyzed. A uniform stress a x of 1 
psi is applied along the boundary at x = ± 1 in. The material properties are 

E = 30 x 10 6 psi and v = 0.336. The finite element mesh pattern used has 10 

(^-direction) by 9 (r-direction) by 6 (z-direction) mesh divisions with a 
total of 2310 degrees of freedom. The resulting circumferential stresses a <j, 
at z = 0.025 in. are plotted in Figure 17. Also shown in the figure are the 
corresponding 2-D results. It is seen that these two solutions are in good 
agreement with the maximum error being 3% at 90°. 

The second study analyzes a (90/0) s boron/epoxy laminate with a circular 
hole in the center. The plate dimensions and the finite element mesh geometry 
are identical with those shown for the isotropic plate. In order to compare 
the solution material properties given in [11] are used. A uniaxial average 
stress cr x of 1 psi is applied at a remote boundary. The resulting in-plane 
tangential stresses in the 0° ply are shown in Fig. 18 along with those 
obtained in [11-13]. The stresses in the 90° ply are also shown in Fig. 19. 
Good agreement between the present and the other solutions is observed. Fig. 
20 compares the normal stress cr z in the mid-plane from various solutions 
[11,13,14]. On a relative basis, there are more discrepancies in the o z 
stress distribution than the ° $ stress among these solutions. Higher stresses 
are obtained in [11] by the use of special hybrid elements around a circular 
hole. On the other hand, the solution in [13] uses a boundary layer method 
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and is significantly different from others [11,12] near cf> = 90°. In any case, 
some experimental studies remain to be done in order to verify the accuracy of 
the 3-D analytical results. 

B. Viscoelastic Response of Cross-Ply Laminates 

The time-dependent 3-D stresses around a circular hole in graphite/epoxy 

(0/90) s and ( 90/0 ) s laminates are analyzed. A uniform displacement of u x = 

0.005 in is prescribed along the remote boundary x = ± 1.0 in. Due to 

symmetry, only one-eighth of the laminate is considered in the analysis. As 

shown in Fig. 21, this model contains 440 3-D solid elements with 2010 degrees 

of freedom. The material properties given in Table 1 are used. The 

viscoelastic response in the 140°F environment is investigated over a time 

7 

period up to 1 year (3.15 x 10 seconds). The distributions of interlaminar 
normal and shear stresses around the hole edge are obtained as a function of 
time. 


Fig. 22 shows the normal stress o z around a hole in the midplane (z = 0) 
of the (0/90) s laminate. The maximum a z which occurs approximately at <f> = 36° 
relaxes by 28% after one year. For the ( 90/0 ) s laminate the a z distributions 
are completely different as shown in Fig. 23; the maximum stress in the 
(90/0) s layup occurs at <j> = 90°. Furthermore, the magnitude of o z in the 
(90/0) s is one order higher than that of the (0/90) s laminate. Thus, the 
normal a z stress is strongly dependent on the laminate stacking sequence. 

The interlaminar normal stress at the interface (z/h = 1.0) of the 
(0/90) s laminate is shown in Fig. 24. For an applied strain of 0.005 in/in, 
the maximum elastic stress a 2 at the 0/90 interface is about 3700 psi which is 
close to 75% of the static ultimate strength of a matrix material. As a 
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result, del ami nation at this interface is expected to occur as the loading is 
increased. The o z distribution at the interface in the ( 90/0 ) s layup is shown 
in Fig. 25. In the (90/0) s layup the interface stress a z is also large but it 
is smaller than the corresponding midplane stress. In both (0/90) s and 
( 90/0 ) s laminates, the interlaminar stress a z is in tension and can cause 
delaminations between plies. 

The interlaminar shear stresses T z< ^ at the interface in both (0/90)s and 
(90/0)s laminates are shown in Figs. 26 and 27. Approximately opposite 
distributions are observed for the interlaminar shear stress in these two 
1 ami nates. 

The effects of time and temperature on the stress relaxation in a 
cross-ply laminate is illustrated in Fig. 28. The maximum interlaminar normal 
stress at 4> = 90° in the {0/90 ) s laminate is plotted as a function of time at 
70°F and 140°F. In this figure, the stress o z has been normalized with 
respect to its elastic solution at time t = 0. Similar curves are shown in 
Fig. 29 for the (90/0) laminate. After 1 year the maximum stress relaxes 48% 
in the (0/90) s laminate and 37% in the (90/0) s laminate. The normalized 
interlaminar shear stress at 4» = 54° is shown in Fig. 30. Less relaxation is 
observed for the t Z( j, stress than for the CT Z stress. 

C. (45/-45) s Laminates 

In (45/-45) s laminates the conditions of symmetry used for the (0/90)s 
layup are no longer valid. However, other symmetric and antisymmetric 
conditions in displacements and stresses can be utilized so that only a 
quarter of the plate is sufficient for the analysis [15], Fig. 31 shows a 
typical finite element model for the (45/-45) s laminate which has 512 solid 
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elements with 2310 degrees of freedom. The boundary conditions and material 
properties are the same as those used for a (0/90) s laminate. 

The normal stress o z at the midplane is shown in Fig. 32 for T = 140°F. 
It is noted that the a z is distributed between a tensile stress of 4750 psi 
and a compressive stress of 4200 psi with the maximum stress occuring at <p = 
78.8°. At <|> = 33.8° and -56.2° the stress values appear to be independent of 
time. Fig. 33 illustrates the interlaminar shear stress distribution at the 
45/-45 interface. For the applied strain of 0.005 in/in the maximum value of 

T z< p at 4> = 90° is 8840 psi which is large enough to initiate delamination. 

The inplane circumferential stress in the -45° ply is shown in Fig. 34. As 
can be seen in these figures the (45/-45) s layup exhibits much stronger 
viscoelastic effect than a cross-ply laminate. 

D. Quasi -isotropic Laminates 

Three-dimensional analyses are also performed for a quasi-isotropic 
(45/0/-45/90) s laminate with a circular hole at T = 140°F. A quarter of the 
laminate including 1025 elements and 4158 degrees of freedom is used in the 
analysis. The finite element mesh pattern used is similar to that shown in 
Fig. 31 except there are 8 subdivisions in the thickness (z) directions. An 
average a x -stress of 1 psi is applied to the laminate edges x = ± 1 in. This 

condition results in a a x of 2.296 psi, 0.763 psi, and 0.178 psi in the 0, 45 

(or -45) and 90 plys, respectively. These ply stresses are used as boundary 
conditions for the problem. 

The circumferential strains around the hole (r/a = 1.04) in each ply 
are shown in Fig. 35. The ply strain is also compared with the 2-D results 
from Section 5.2. The strain is fairly uniform throughout the laminate 
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thickness. The time-dependence of the £ <j> strain in the 90° layer is shown in 
Fig. 36. A maximum of 14.6% increase in £ <j> is obtained after one year. It is 
interesting to note that the £ <j, distribution in each layer is almost symmetric 
about the centerline <t> = 0. This suggests that additional symmetry conditions 
can be imposed in the finite element model, i.e. only one eighth of the 
laminate is sufficient for the calculation of e<j, in a quasi -isotropic 
1 aminate. 

The transverse strain £ z distributions are plotted in Fig. 37 for the 
90-ply (z/h = 0.75) and in Fig. 38 for the 0-ply (z/h = 2.75). These 
distributions are quite different from the £ ( j ) results; the transverse strain 
£ z is not symmetric with respect to 4> = 0. Furthermore, the amount of strain 
increase in the 0-ply is relatively large. An explanation for this behavior 
is that during creep the 45 and 90 plies lose stiffness since the material 
properties degrade most in these layers. This yields much higher strains in 
the 45 and 90 plies than in the 0-ply. However, due to the compatibility 
condition between adjacent layers additional interlaminar shear stresses are 
built up and are added to the 0-ply. As a result, the actual loads increase 
in the 0-ply during creep resulting in a much larger strain than would be 
obtained for a unidirectional composite. This load transferring mechanism is 
illustrated in Fig. 39 in which the in-plane stress a x near the remote 
boundary x = 1.0 is shown at two different times. The a x stress increases 10% 
in the 0-ply but decreases 15.6% in the 45-ply and 32% in the 90-ply after one 
year. 


Figure 40 depicts the interlaminar shear strain y Z( p in the -45 ply around 
a circular hole. This distribution is a mirror image of that in the +45 ply. 
Over a period of one year the rate of increase in shear strain is 60% at 
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<j> = +62° in the -45° ply. The high creep strain rate is due to the fact that 
the GY70/339 material used in the analysis has strong time-dependent 
properties. 

E. Matrix Dominated Laminates 

The last example analyzed is a ( 90/-45/90-45 ) s laminate whose properties 
are dominated by the matrix material. The boundary conditions applied are a x 
= 0.389 psi in the 90-ply and a x = 1.621 psi in both 45 and -45 plies. These 
conditions yield an average stress a x of 1 psi across the laminate thickness. 

Results of strains e<j, in each ply of the ( 90/-45/90-45 ) s laminate at T = 
140°F are compared with the corresponding 2-D solution in Fig. 41. The strain 
component e z , and y<j, z in the 90 layer (z/h = 3.75) are plotted in Figs. 

42, 43, and 44, respectively. As expected, much higher creep rate is observed 
in a matrix dominated layup than in a quasi -i sotropic laminate. Specifically, 
the maximum strain increases 50% in the (90/-45/90/45)s layup as compared 
to 14.6% in the (45/0/-45/90)s layup over a time period of one year. As 
before, such a large creep rate is due to the material properties assumed in 
the analysis. 
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6. CONCLUSION 


The time dependent behavior of composite materials has been described by 
an anisotropic thermo-viscoelastic constitutive model. A numerical procedure 
has been developed for the solution of time-dependent stresses and strains in 
composite laminates containing geometric discontinuities and complicated 
boundary conditions. Using this procedure, the stress and strain 
distributions around a circular hole in graphite/epoxy composites have been 
obtained as a function of time for both mechanical and thermal loads. The 
effects of layup orientation and load spectrum on deformation histories have 
been demonstrated. The results show that the present method gives the 
accurate and efficient numerical solutions of complex anisotropic 
thermoviscoelastic boundary value problems. 
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APPENDIX A 


1. For the two dimensional analysis, 
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and all other n-j j r = 0 (i,j = 1,2, ... 6, and r = 1,2, ... 9). 
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u 

x i 


0.06698253 



0.0729459 

8.174141919E+15 


0.0696426 

4.976486103E+14 


0.150514 

1 . 477467 149E+ 13 


0.148508 

4.761315266E+11 

" 

0.146757 

1.799163029E+10 


0.102892 

5.253922053E+08 


0.114155 

1.846670914E+07 


0.071036 

5.288067476E+05 


0.0484272 

1.494783951E+04 


0.00813977 

5.516602214E+02 
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Table 3 


Shift Factors for Various Temperatures 


Temperature (°F) 

Shift Factor 

i 

75 

1.0 

104 

8.9125 

122 

7.9433E+1 

140 

1.5849E+3 

160 

6.6069E+4 

212 

6.3096E+11 

250 

1.0000E+18 


.Table 4 Elastic Stress Concentration Factors (S.C.F.) 



Coarse 

Mesh 

Intermediate 

Mesh 


Ref. [10]* 

S.C.F. 

3.2470 

3.2625 

3.2514 

3.2292 

No. of Elements 

110 

297 

570 


No. of DOF's 

268 

670 

1246 



*The solution has been corrected for the finite width effect. 
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Fig. 1 Coordinate systems for a unidirectional composite 



Fig. 2 Discretization of temperature history 


28 


Fig. 3a Finite element mesh pattern for a unnotched laminate (2-D) 



Fig. 3b Finite element mesh pattern for a notched laminate (2-D) 


29 






Fig. 4 Normalized relaxation function 
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Fig. 5 In-plane strains e x in a (0/45/90/-45) s laminate 
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Fig. 16 Average circumferential stresses around a hole in (45/-45) s laminate 
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Fig. 17 Comparison of the 3-D and 2-D solutions 
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Fig. 18 Elastic a<j) stresses in the 0° layer of a (90/0) s boron/epoxy laminate 
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Fig. 19 Elastic CT(j> stresses in the 90° layer of a (90/0) s boron/epoxy laminate 
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Fig. 21 Finite element mesh pattern used for cross-ply laminates 
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Fig. 22 Interlaminar normal stress Cz at the mid-plane in a (0/90) s laminate 
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Fig. 25 Interlaminar normal stress a z at the interface in a (90/0) s laminate 
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Fig, 26 Interlaminar shear stress x z <j> at the interface in a (0/90) s laminate 
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Fig. 27 Interlaminar shear stress T Z( j) at the interface in a (90/0) s laminate 
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Fig. 28 History of the interlaminar normal stress at the interface of a (0/90) s laminate 
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Fig. 29 History of the interlaminar normal stress at the interface of a (90/0) s laminate 
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Fig. 30 History of the interlaminar shear stress at the interface of a (0/90) s laminate 
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Fig. 33 Interlaminar shear stress % z § at the interface in a (45/-45) s laminate 
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Fig. 36 Circumferential strain e<j> in the 90° ply of a (45/0/-45/90) s laminate 
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Fig. 39 Through-the-thickness stress distributions ( x/L = 0.9, y/W = 0.125 ) 
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Fig. 41 Comparison of elastic ply strains e<j> with the 2-D solution for a (90/-45/90/45) s laminate 
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Fig. 42 Circumferential strain e<s in the 90° ply of a (90/-45/90/45) s laminate 
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Fig. 43 Transverse strain e z in the 90° ply of a (90/-45/90/45) s laminate 
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of a (90/-45/90/45) s laminate 



